#include <cstdlib>
#include <climits>
#include <cstdint>
#include <cassert>

#include <iostream>
#include <string>
#include <vector>
#include <functional>
#include <unordered_map>
#include <exception>
#include <stdexcept>
#include <algorithm>

using std::vector;
using std::unordered_map;
using std::string;

enum class Nucleotide {
	A, C, G, TU
};

struct ScoringValues {
	vector<vector<int> > matrix;
	int opening_penalty;
	int extension_penalty;
	int threshold;
};

struct Block {
	size_t start;
	size_t end;

	bool operator==(const Block& b) const {
		return start == b.start && end == b.end;
	}
};

namespace std {
template<>
struct hash<Block> {
	size_t operator()(const Block& b) const {
		return b.start;
	}
};
}

class SplicedAlignment {
public:
	SplicedAlignment(const vector<Nucleotide>& dna,
			const vector<Nucleotide>& rna, const vector<Block>& blocks,
			const ScoringValues& scoring_vals) :
			dna(dna), rna(rna), blocks(blocks), scoring_vals(scoring_vals) {}

	void Compute() {
		auto rit = blocks.rbegin();
		while (rit != blocks.rend()) {
			const Block& b = *rit;
			Compute(b.end, rna.size() - 1, rit, false);
			rit++;
		}
	}

	const vector<Block> GetResult() const {
		vector<Block> res;
		const Block *block = FindOptimalEndBlock();
		assert(block != nullptr);
		res.push_back(*block);
		// traverse linked list of blocks to the start
		int j = rna.size() - 1;
		while (j >= 0) {
			try {
				// jumps[j] may not exist
				const Block *b = jumps.at(j).at(*block);
				if (b != nullptr) {
					res.push_back(*b);
					block = b;
				}
			} catch (std::out_of_range& e) {
				// continue
			}
			j--;
		}
		std::reverse(res.begin(), res.end());
		return res;
	}

private:
	const vector<Nucleotide>& dna;
	const vector<Nucleotide>& rna;
	const vector<Block>& blocks;
	const ScoringValues& scoring_vals;
	vector<vector<unordered_map<Block, int> > > pts;
	// indexed by j and block to the previous block, nullptr if we don't jump further
	vector<unordered_map<Block, const Block*> > jumps;


	// open_gap specifies whether at this position a gap is already open (i.e. we have already
	// paid the gap opening penalty or not)
	int Compute(size_t i, size_t j, vector<Block>::const_reverse_iterator it, bool open_gap) {
		// exit condition
		if (i == 0 || j == 0) {
			int val = std::max(i, j);
			int start_penalty = open_gap ? scoring_vals.extension_penalty : scoring_vals.opening_penalty;
			return (val - 1) * scoring_vals.extension_penalty + start_penalty;
		}
		const Block& b = *it;
		try {
			return pts.at(i).at(j).at(b);
		} catch (std::out_of_range& e) {
			// value does not exist yet so compute it now
		}
		// if gap is already open we only pay extension_penalty instead of opening_penalty
		int gap_penalty = open_gap ? scoring_vals.extension_penalty : scoring_vals.opening_penalty;
		int match_value = scoring_vals.matrix[(int)dna[i]][(int)rna[j]];
		if (i == b.start) {
			// We possibly have to find a new block now to continue the alignment
			int max_val = Compute(i, j - 1, it, true) + gap_penalty;
			const Block *next_block = nullptr;
			// blocks are sorted by start position, so all blocks before current one may
			// precede it, but none after it
			while (it != blocks.rend()) {
				const Block& candidate = *it;
				if (candidate.end <= b.start && ValidCandidate(candidate, j, open_gap)) {
					int val = Compute(candidate.end, j, it, true) - gap_penalty;
					if (val >= max_val) {
						max_val = val;
						next_block = &candidate;
					}
					val = Compute(candidate.end, j - 1, it, false) + match_value;
					if (val >= max_val) {
						max_val = val;
						next_block = &candidate;
					}
				}
				it++;
			}
			AssignToPts(i, j, b, max_val);
			AssignToJmps(j, b, next_block);
		} else {
			AssignToPts(i, j, b, std::max(
					std::max(
							Compute(i - 1, j, it, true) - gap_penalty,
							Compute(i, j - 1, it, true) - gap_penalty),
							Compute(i - 1, j - 1, it, false) + match_value));
			AssignToJmps(j, b, nullptr);
		}
		return pts[i][j][b];
	}

	// Checks whether the given block is a valid candidate at the gene position j
	// with the given threshold, i.e. that its local score is at least the threshold
	bool ValidCandidate(const Block& block, int j, bool open_gap) const {
		return GetAlignmentScore(block.start, block.end, j, open_gap) >= scoring_vals.threshold;
	}

	// Compare alignment until min_i is reached
	int GetAlignmentScore(int min_i, int i, int j, bool open_gap) const {
		// if gap is already open we only pay extension_penalty instead of opening_penalty
		int gap_penalty = open_gap ? scoring_vals.extension_penalty : scoring_vals.opening_penalty;
		int match_value = scoring_vals.matrix[(int)dna[i]][(int)rna[j]];
		if (i == min_i) {
			return std::max(gap_penalty, match_value);
		}
		return std::max(GetAlignmentScore(min_i, i - 1, j - 1, false) + match_value,
				std::max(GetAlignmentScore(min_i, i - 1, j, true) - gap_penalty,
						GetAlignmentScore(min_i, i, j - 1, true) - gap_penalty));
	}

	void AssignToPts(size_t i, size_t j, const Block& b, int value) {
		if (pts.size() <= i) {
			pts.resize(i + 1);
		}
		if (pts[i].size() <= j) {
			pts[i].resize(j + 1);
		}
		pts[i][j][b] = value;
	}

	void AssignToJmps(size_t j, const Block& b, const Block* target) {
		if (jumps.size() <= j) {
			jumps.resize(j + 1);
		}
		jumps[j][b] = target;
	}

	const Block* FindOptimalEndBlock() const {
		const Block* end_block = nullptr;
		int max_val = INT_MIN;
		std::for_each(blocks.begin(), blocks.end(), [&](const Block& block) {
			// value has to exist here for endblock!
			int val = pts[block.end][rna.size() - 1].at(block);
			if (val > max_val) {
				end_block = &block;
				max_val = val;
			}
		});
		return end_block;
	}

};

vector<Nucleotide> ParseString(string data) {
	vector<Nucleotide> res;
	std::for_each(data.begin(), data.end(), [&res](const char ch) {
		switch(tolower(ch)) {
		case 'a':
			res.push_back(Nucleotide::A);
			break;
		case 'c':
			res.push_back(Nucleotide::C);
			break;
		case 'g':
			res.push_back(Nucleotide::G);
			break;
		case 'u':
		case 't':
			res.push_back(Nucleotide::TU);
			break;
		default:
			throw std::runtime_error("Illegal character in string");
		}
	});
	return res;
}

int main(int argc, char **argv) {
	vector<Nucleotide> dna = ParseString("ATCG");
	vector<Nucleotide> rna = ParseString("AAAA");
	vector<Block> blocks {
		{0,3}
	};
	// Scoring matrix
	vector<vector<int> > matrix {
		{1, 0, 0, 0},
		{0, 1, 0, 0},
		{0, 0, 1, 0},
		{0, 0, 0, 1},
	};
	ScoringValues scoring_values {
		matrix,
		1, // opening_penalty
		1, // extension_penalty
		0  // threshold
	};
	auto code = SplicedAlignment(dna, rna, blocks, scoring_values);
	code.Compute();
	auto result = code.GetResult();
	std::for_each(result.begin(), result.end(), [](const Block& block) {
		std::cout << "(" << block.start << "," << block.end << ")" << std::endl;
	});
	return EXIT_SUCCESS;
}

